#include <src/math/norms.h>
#include <src\algorithm\Evolution_Algorithm\TimeStepAlgorithm.h>

#if !defined(_CrankNicolson_)
#define _CrankNicolson_


template <class T>
struct CrankNicolson_Interface
{
	typedef CrankNicolson_Interface								type;
	CrankNicolson_Interface() {};

	virtual OS_double CFL() = 0;
	virtual bool criteria() = 0;
	virtual void gstar() = 0;
	virtual void vstar() = 0;

};



/*! \brief Crank Nicolson scheme
CrankNicolson generic schemes solved by a fixed point algorithm.
Given the initial condition S^{t), the time increment \Delta, and a non linear operator Op(S),
the CrankNicolson scheme is the solution of the following non linear equation 

S^{n+1} = S^{t) + \Delta_t Op(t^{1/2} ,S^{n+1/2) ), where S^{n+1/2) = \frac {S^{t)+S^{n+1}}}{2} and t^{1/2} = t + \Delta_t /2.

the fixed point algorithm try to solve the previous equation by the following algorithm

S^{0,n+1} : = S^{t); 
k>=1, loop
S^{k,n+1/2) := \frac {S^{t)+S^{k-1,n}{2} 
S^{k,n+1} : = S^{t) + \Delta_t Op(t^{1/2} ,S^{k,n+1/2) ) ;

until a convergence holds (Predictor), default being the \ell^2 norm

|S^{k+1,n+1}-S^{k,n+1}|_{\ell^2}^2 < \epsilon, where |S|_{\ell^2}^2 = \frac{1}{2N} \sum_{i=1,..N} |S|_i ^2

This algorithm always converges for \Delta small enough, under reasonnable conditions over the Operator Op.
However, much care has to be taken for the choice of the CFL condition, and we embed a simple time adaptive technique :
We verify \Delta \le cfl at the beginning of the algorithm.
If the algorithm seems to diverge (more than 10 iterations), then we try to divide the step time \Delta by 2 and reset the algorithm.
If the algorithm resumes to the Forward Euler Algorithm (only one iteration), then the step time is multiplied by two.

*/


template <	class TimeTypeOperator >
class CrankNicolson : public TimeStepAlgorithm<TimeTypeOperator>
	,public CrankNicolson_Interface<TimeTypeOperator>
{
public:
	typedef typename TimeStepAlgorithm<time_operator_type>			Base;
	typedef typename CrankNicolson<time_operator_type>				algorithm_type;
	typedef typename TimeTypeOperator::argument_type				argument_type;
	typedef typename TimeTypeOperator::argument_base_type				argument_base_type;
	virtual ~CrankNicolson(){};
	CrankNicolson() : cfl_adjust_(1.), n_(0) {};
	virtual void gstar() {
		Op_iterate = Op_(t_12,S_k12);
	};
	virtual void vstar() {
		S_k12   = (*get_state_variable()+iterate)/2.;
	};

	result_type operator()(const time_type & from, time_type & to, argument_type S_k) {

		nb_iterations_ = 0.;
		delta_t_	= to-from ;
		iterate		= S_k;
		S_k12		= S_k;
		t_12        = from;
		Op_iterate  = Op_(t_12,S_k12);
		delta_t_ =  std::min(delta_t_, cfl_adjust_*CFL());
//		OS_DYNAMIC_CHECK(delta_t_ > THRESHOLD,"CrankNicolson : the time step is too low. Check the CFL function.");
		vstar();
		gstar();
		iterate    = S_k+delta_t_*Op_iterate;
		t_12		= from + delta_t_ /2. ;
		bool Convergent = false;
		do {
start:				
				nb_iterations_++;

				Op_S_k1 = Op_iterate;
				S_k1    = iterate;
				vstar();
				gstar();
				iterate = S_k + delta_t_*Op_iterate;

				if ( criteria() ) {
					break;
				}
				if (nb_iterations_ > 10 ) {
					adjust_cfl(from,S_k);
				}
		} while (true);
		if (nb_iterations_ < 3) {
//From one side, we should never overpass the cfl condition, that ensures stability.
//			if (cfl_adjust_ < 1)		
				cfl_adjust_ *= 2.; // multiply CFL constant
//From other side, it may be really hard to determine a correct CFL function. 
// May be a external parameter defaulted to false is a good choice.

		}
end:
		to = from + delta_t_;
		++n_;
		return iterate;
	}
	static std::string get_name(){return " Crank Nicolson \delta_t u = ";};
protected :
	virtual void adjust_cfl(const time_type & from, argument_type S_k){
		OS_DYNAMIC_CHECK(cfl_adjust_ > THRESHOLD,"CrankNicolson : problem in Crank Nicolson scheme");
		cfl_adjust_ /= 2.; // divide CFL constant
		delta_t_ /= 2.; 
		t_12		= from + delta_t_ / 2.;
		iterate = *get_state_variable();
		vstar();
		gstar();
		iterate =  (S_k+delta_t_* Op_iterate);
		nb_iterations_ = 0;
	}
	indice_type n_;
	OS_double cfl_adjust_;
	time_type t_12, delta_t_;
	argument_base_type S_k12, iterate, Op_iterate, S_k1, Op_S_k1;
	std::size_t nb_iterations_;
};


#endif
